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(N ■ ABSTRACT 
D 

pLn ', The flux anomalies in four-image gravitational lenses can be interpreted as evidence 

for the dark matter substructure predicted by cold dark matter (CDM) halo models. 
In principle, these flux anomalies could arise from alternate sources such as absorption, 
^ ■ scattering or scintillation by the interstellar medium (ISM) of the lens galaxy, problems 

. in the ellipsoidal macro models used to fit lens systems, or stellar microlensing. We 

, apply several tests to the data that appear to rule out these alternate explanations. 

[ First, the radio flux anomalies show no significant dependence on wavelength, as 

^T) I would be expected for almost any propagation effect in the ISM or microlensing by the 

[ stars. Second, the flux anomaly distributions show the characteristic demagnifications 

of the brightest saddle point relative to the other images expected for low optical 
depth substructure, which cannot be mimicked by either the ISM or problems in the 
macro models. Microlensing by stars also cannot reproduce the suppression of the 
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c/3 . bright saddle points if the radio source sizes are consistent with the Compton limit 

^ ' for their angular sizes. Third, while it is possible to change the smooth lens models 

to fit the flux anomalies in some systems, we can rule out the necessary changes in 
^ ' all systems where we have additional lens constraints to check the models. Moreover, 

' the parameters of these models are inconsistent with our present observations and 

expectations for the structure of galaxies. We conclude that low-mass halos remain the 
best explanation of the phenomenon. 



Subject headings: cosmology: gravitational lensing 



1. Introduction 

In Dalai & Kochanek (2002, DK02 hereafter) we demonstrated that the incidence of anomalous 
flux ratios in gravitational lenses was consistent with the expected mass fraction of satellites found 
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in cold dark matter (CDM) halo simulations (e.g. Kauffmann 1993, Moore et al. 1999, Klypin et 

al. 1999, Bode, Ostriker & Turok 2001, Zentner & Bullock 2002). Our quantitative estimate for the 
substructure mass fraction was part of an outburst of interest in substructure and anomalous flux 
ratios of specific lenses (e.g. Mao & Schneider 1998, Keeton 2001, Bradac et al. 2002, Chiba 2002, 
Metcalf 2002), in CDM models (e.g. Metcalf & Madau 2001, Chiba 2002, Zentner & Bullock 2002), 
and the difficulty in explaining them based on the properties of the primary lens galaxy (e.g. 
Metcalf & Zhao 2001, Keeton, Gaudi & Fetters 2002, Evans & Witt 2002, Moeller, Hewett & 
Blain 2002, Quadri, Moeller &: Natarajan 2002). In DK02 we argued that alternate explanations 
to substructure for the anomalous flux ratios were unlikely, but we did not explore the alternate 
possibilities in detail or develop tests to distinguish substructure from the alternatives. We do so 
in the current paper. 

The alternatives to substructure fall into three broad categories: propagation effects in the 
interstellar medium (ISM) of the lens galaxy, problems in the "macro" models for the gravitational 
potential of the lens galaxy, and confusing "microlensing" by the stars in the lens galaxy with the 
effects of more massive satellites. In this paper we will consider all three of these possibilities. 
In §2 we study the wavelength dependence of the flux anomalies and develop a simple statistical 
test to demonstrate that the anomalous flux ratio problem must be due to gravity. In §3 we show 
that deviations from elliptical shapes appear unable to account for flux anomalies, and that errors 
in the model for the potential of the primary lens have statistical properties differing from those 
created by substructure and observed in the data. In §4 we show that it is difficult to explain 
the anomalous flux ratios of the radio lenses using microlensing. We summarize our results and 
outline further tests in §5. 

2. Ruling out the ISM 

The ISM can affect flux ratios either through absorption or scattering. In optical/IR 
observations, many gravitational lenses show wavelength-dependent flux ratios consistent with 
the effects of dust extinction (e.g. Falco et al. 1998). In the radio, some lenses show the effects 
of scatter broadening by electrons at low frequencies (e.g. B1933+503, Marlow et al. 1999a), 
although in none of these systems is there evidence for a net change in the radio flux. In this 
section wc examine whether the radio lenses show any frequency dependent changes in their flux 
ratios that could be a signature for the ISM modifying the observed fluxes, and develop a simple, 
non-parametric test to distinguish the effects of the ISM from the effects of gravity. 

2.1. The Frequency Dependence of Anomalous Flux Ratios 

Almost all mechanisms by which the ISM can modify flux ratios should show frequency- 
dependent effects, since most scattering processes relevant to the fluxes of radio images become 
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weaker at higher frequencies. For example, weak scintillation causes flux perturbations with a 
scaling (assuming Kolmogorov turbulence) oc i/"^''/^^ (Narayan 1992) and the optical depth for 
free-free absorption scales roughly as oc z/~^-^ (e.g. Mczgcr & Henderson 1967). Thus, one way 
to constrain the ISM as an explanation of the anomalous flux ratios is to examine the allowed 
frequency dependence of the effect. If the source has an intrinsic spectrum f^^g as a function of 
frequency u, then in the absence of any perturbations from the ISM, the images have spectrum 
fi^^i = |Mj|/,y,s given the (signed) image magnification Mj. If the ISM modifies the fluxes through 
a frequency-dependent optical depth Ti^^i, then the image fluxes become fi^^i = |Mj|/i^^s exp(— T,/^j). 

We examined the optical depth function using the four (or more) image radio lenses 
with fluxes measured at 5 GHz and either 8 GHz or 15 GHz (MG0414+0534, Katz, Moore 
k Hewitt 1997; B0712+472, Jackson et al. 1998; B1359+154, Myers et al. 1999; B1422+231, 
Patnaik & Narasimha 2001; B1555+375, Marlow et al. 1999b; B1608+656, Fassnacht private 
communication; B1933-F503, Sykes et al. 1998; and B2045-I-265, Fassnacht et al. 1999). We used 
the same models and magniflcation estimates as were used to estimate the substructure fraction 
with a 5% lower bound on the flux errors. The simplest way to illustrate the problem is to assume 
that one image is strongly affected by the ISM while the others suffer from little or no absorption 
or scattering. For each image in a lens, we used the fluxes of the other three to estimate the 
intrinsic spectrum of the source. From this estimate of the intrinsic spectrum we computed the 
optical depth needed to reproduce the observed flux of the remaining image by computing the 
necessary optical depth at 5 GHz, rs, and its spectral index, a, under the assumption that r oc z/"^ 
between 5 GHz and the available higher frequency. We estimated the errors solely from the 
published measurement errors, so they include no uncertainties arising from the lens model or any 
time variability between the measurement epochs for the different frequencies. 

The results are shown in Fig. 1. We have coded the images by their parities (minima versus 
saddle points) and relative fluxes (brightest versus faintest for each parity) because we will 
show below that the distributions of the model flux residuals depend strongly on these image 
identifications. This is seen in Fig. 1 as the concentration of the bright saddle-point images 
in the direction of higher optical depth compared to the other images. We have dropped four 
(of 32) images where the estimated optical depth changes sign between the higher and lower 
frequencies and it is impossible to define a spectral index. Except at zero optical depth, where the 
slope estimates become unstable and have larger uncertainties, the lenses broadly require a ~ 0. 
Since most radio propagation effects have strong wavelength scalings, the ISM seems an unlikely 
explanation. 



2.2. A Statistical Test for Substructure 

While the spectral indices required to explain the anomalous flux ratios are peculiar for 
standard scattering or absorption processes, we would prefer to have a test which can distinguish 
any ISM effect from the gravitational effects of substructure. One approach is to note that the 
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Fig. 1. — Estimates of the ISM properties needed to explain the anomalous flux ratios. We show 
an estimate of optical depth at 5 GHz, rs, and the spectral index a, where r oc i/", needed to 
bring the observed flux of each image into agreement with the average source properties predicted 

by the other three images. The images are coded by their parities and relative fluxes with squares 
used for minima (positive parity) and triangles for saddle points (negative parity). Filled symbols 
are used for the brightest image of each parity, and open symbols for the faintest image. Note the 
concentration of the brightest saddle points (minima) towards positive (negative) optical depth. 
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physical properties of the ISM are a local property of the lens, and should have no knowledge of 

the global properties of the lens geometry. Radial gradients in the properties of the ISM should 
matter little because in each lens the images of interest effectively lie at the same radius and 
there is no simple angular correlation between the positions of the images and the major axis of 
the galaxy. Moreover, creating the anomalous flux ratios depends on the presence of clumped 
ISM components rather than smoothly distributed components both to create physical conditions 
extreme enough to have an effect and to differentially affect the lensed images. The only exception 
to this rule arises for effects such as scintillation or scatter broadening where the effect diminishes 
as the source size becomes larger - for these cases, the ISM should preferentially perturb the least 
magnified images because they are the most compact. 

The effects of gravity, however, are not determined purely by local conditions because the 
image positions and fluxes are determined by the competition between the local gravitational 
potential and the geometric time delay. The magnification tensor, whose eigenvalues determine the 
image parities, depends on the projected surface density, the projected tidal shear in the gravity, 
and the redshifts of the lens and the source. As a result, we can attach a global identification 
to each image based on its parity (maximum, minimum, saddle point of the time delay surface) 
as well as its magnification with which the ISM should show little (magnification) or no (parity) 
correlation. Image parities and magnification orderings (but not the precise magnifications) are 
generic predictions of successful lens models. In particular, for a four image lens, the images 
alternate parities in the sequence minimum-saddle point-minimum-saddle point as we go around 
the lens. We can also identify the most and least magnified image of the two minima or saddle 
points from the geometry of the lens. Unlike the ISM, the effects of substructure depend on the 
image parity and magnification. The fluxes of highly magnified images arc more unstable to small 
gravitational perturbations than the fluxes of less magnified images (as originally emphasized by 
Mao & Schneider 1998). Most importantly, the magnification perturbations created by low optical 
depth substructure have a very strong dependence on the image parity because the perturbations 
to the fluxes of the saddle points arc skewed in the direction of demagnification (Schcchtcr & 
Wambsganss 2002, Keeton 2002). This leads to a simple non-parametric test for distinguishing 
the effects of the ISM from the effects of gravity - if the flux anomalies depend on parity and 
magnification they cannot be due to the ISM. 

We first demonstrate the effects of substructure on the distribution of flux residuals. We took 
the same sample of 7 lenses we used in DK02 (MG0414+0534, Hewitt et al. 1992; B0712+472, 
Jackson et al. 1998; PG1115+080, Weymann et al. 1980; B1422+231, Patnaik et al. 1992; 
B1608+656, Fassnacht et al. 1996; B1933+503, Sykes et al. 1998; and B2045+265, Fassnacht et 
al. 1999). The data were fit using standard macro models (singular isothermal ellipsoids (SIE) 
plus external shear) to generate an initial model. We then added a substructure mass fraction 
fsat near the images, modeling the substructure as tidally truncated isothermal spheres with 
critical radius b = O'.'OOl (corresponding to approximately 10® M©). With the substructure added 
we reidentified the lensed images and their fluxes and added astrometry (0'.'003 rms) and flux 
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Fig. 2. — Cumulative distributions of model flux residuals, log{fobs/fmod), expected from Monte 
Carlo simulations of substructure with mass fraction fsat between 0.20 (top) and 0.002 (bottom). 
The solid (dashed) lines are for minima (saddle points), with squares (no squares) for the 
distribution corresponding to the most (least) magnified image. Note the steady shift of the 
brightest saddle point to fainter flux residuals as the satellite fraction increases. For very large 
satellite fractions the distributions would become more similar. The Monte Carlo samples are 
based on the same lens sample as was used in DK02 assuming 10% flux measurement errors. 
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Fig. 3. — Cumulative distributions of model flux residuals, log{fobs/fmod), in the real data, assuming 
constant fractional flux errors for each image. The solid (dashed) lines are for minima (saddle 
points) , with squares (no squares) for the distribution corresponding to the most (least) magnified 
image. Prom top to bottom the distributions are shown for samples of 8 radio, 10 optical or 15 total 
four-image lenses. If the flux residuals are created by propagation effects we would not expect the 
distributions to depend on the image parity or magnification, while if they are due to low optical 
depth substructure we would expect the distribution for the brightest saddle points to be shifted 
to lower observed fluxes. Note the similarity between the observed curves and the Monte Carlo 
simulations for substructure fractions of order ~ 5%. 
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(10% rms) measurement errors to generate a model data set with substructure. This model data 
was then refit using the standard macro models, including no weight on the fit to the fluxes 
of the images. We then compare the observed image fluxes, fobs,i to the model image fluxes, 
fmod,i = l-^il/src) found given the magnification Mj predicted by the model fitted to the perturbed 
data. We generated 10 such realizations for each of the 7 lenses. Figure 2 shows the distributions 
of the residuals, ^og{ f obs,i/ fmod,i)j for individual image types (saddle point/minimum, 
brightest/faintest) as a function of the substructure mass fraction, fsat- As expected, the residual 
distribution for the brightest saddle point is markedly different from that for all other images 
unless the substructure fraction is so low as to be masked by the random flux errors. 

We now repeat the test for the real data. The same standard macro models were used for 
each system, and we again assigned no weight to fitting the observed image fluxes. For radio lenses 
we used the highest frequency measurements from either the VLA or Merlin, and for the optical 
lenses we used the flux ratios at the longest available wavelength (generally H-band where none 
of the systems shows significant extinction). This gave us a sample of 8 radio systems (adding 
B0128+437, Phillips et al. 2000, and B1555+375, Marlow et al. 1999b to the DK02 sample), 10 
optical systems ^ and 15 joint systems. For our standard model we estimated the source flux using 
constant fractional errors on the image fluxes. Figure 3 shows the distribution for the residuals, 
log(/o6s//mod) foi' all presently available quads containing either compact radio sources or bright 
quasars. Several systems appear in both categories, and for these systems we use the radio fluxes 
for the joint sample. The distributions of the residuals differ for minima and saddles in the sense 
that saddle points tend to be fainter than expected and minima tend to be brighter - just as in 
the simulations of the distributions expected for substructure in Fig. 2. 

We used the Kolmogorov-Smirnov (KS) test to estimate the significance of the differences. If 
we simply compare the distributions of log{fabs/ fmod) ^or minima and saddle points, we find that 
the KS test probabilities for the two distributions to be the same are 2.3%, 28%, 5.5% for the 
radio, optical and joint samples respectively. Theoretically, we expect the most highly magnified 
saddle point to have the most discrepant residual distribution because the high magnification 
makes the flux unstable to perturbations and because the theoretical studies of low optical 
depth substructure predict that it should show the largest differences (as in Fig. 2). The KS 
test probabilities for the most magnified saddle point to have the same residual distribution as 
the other three images is only 0.04%, 5% and 0.3% for the radio, optical and joint samples. In 
each case, the most highly magnified minimum shows the next lowest probability for a residual 
distribution agreeing with the other images (1.9%, 59% and 18% respectively), but the significance 
is lower. In the case of substructure, the magnification perturbation distribution for the minima 



^The systems are the four-image quasar lenses HE0230-2130 (Wisotzki et al. 1999), HE0435-1223 (Wisotzki 
et al. 2002), HS0810+2554 (Reimers et al. 2002), RXJ0911+0551 (Bade et al. 1997), PG1115+080 (Weymann et 
al. 1980), H1413+117 (Magain et al. 1988), Q2237+030 (Huchra et al. 1985), B1422+231 (Patnaik et al. 1992), 
B2045+265 (Fassnacht et al. 1999) and MG0414-I-0534 (Hewitt et al. 1992). The latter 3 systems are also in the 
radio sample. We used the Castles H-band fluxes if available, otherwise the I-band fluxes. 
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is harder to differentiate from the measurement uncertainties which probably dominated the 
residuals for the faintest images. 

We can check these significance estimates by examining trials in which the identifications of 
the images are randomized. We made 10^ trials in which we assigned the image identifications 
(l=brightest minimum, 2=faintcst minimum, 3=brightcst saddle point, 4=faintcst saddle point) 
to the images at random with the restriction that each lens always have the 4 possible image types. 
We then estimated the probability that these random assignments could produce results similar 
to that found for the true assignments. The fraction of the time the observed values of the KS 
statistic were exceeded tracked the KS test probability estimates well. For example, in only 2.3%, 
18% and 4.8% of the trials did the KS statistic for the comparison between the distributions of 
minima and saddle points exceed that of the true distribution. We also examined the distribution 
of the KS statistic if we bootstrap resampled the lens sample. The KS statistics and probabilities 
were typical of any sample of lenses drawn (with replacement) from the available 4-image lenses, 
so the results are not dependent on any particular system or small subset of systems. 

We also examined weighting schemes for estimating the source flux other than assuming 
constant fractional errors for the observed image fluxes. Note, however, that our test is independent 
of any simple rescaling of the source flux. The KS test makes a non-parametric comparison of the 
^og{ fobs / f mod) distributions, and a uniform change in the source flux simply shifts the distribution. 
If we assume constant flux errors rather than constant fractional errors, the probability the saddles 
and minima have the same distribution is again 2.3%. If we set the source flux using the fluxes of 
the faintest minimum and saddle, then the probability that the brightest minimum and saddle 
have the same residual distribution for the radio lenses is 0.02%. 

Finally, we can compare the observed distribution to either those expected for Gaussian 
errors or those we calculated for our Monte Carlo simulations of substructure. For example, while 
the overall residual distributions of the images are well fit by a log-normal distribution of width 
(Tin/ = 0.24, the residual distributions of the images separately arc grossly inconsistent with 
log-normal distributions. The KS test likelihood for the joint distribution of all images matching 
a log-normal distribution is 0.26, while the likelihood for the distributions of the individual types 
matching a single log- normal distribution is 1.5 x 10~^. This is simply another way of quantifying 
the fact that the residual distributions depend on the image parities and flux orderings. If we 
compare the residual distributions by image for the radio lens sample to the distributions predicted 
in the Monte Carlo simulations, we find KS test probabilities for the distributions to be the same 
of 4 X 10^^, 2 X 10^^, 2 X 10"'^, 0.02, 0.3 and 0.05 for the cases with satellite mass fractions of 

= 0.002, 0.01, 0.02, 0.05, 0.10 and 0.20 respectively. While we do not, at present, feel this 
is as optimal a way of estimating fsat as the methods we used in Dalai &; Kochanek (2002), it 
emphasizes the consistency of the observed dependence of the residual distribution on parity and 
flux with that expected for substructure. 
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Fig. 4. — Surface density contours for models of PG1115+080 including misaligned 03 and 04 
niultipolcs (thin lines). The model in the left panel is constrained only by the 4 compact images 
(Ai, A2, B and C, filled squares). The model in the right panel is also constrained by the Einstein 
ring image of the host galaxy. The heavy dashed curve is the tangential critical curve of the lens. 
The heavy solid curve is the ring curve of the host galaxy and the thin dashed curve is the best fit 
model. 




-1 -0.5 0.6 1 -1 -0.5 0.5 1 



X (arcsec) x (arcsec) 

Fig. 5. — Surface density contours for models of B1933+503 including misaligned 03 and 04 
multipoles (thin lines). The model in the left panel is constrained only by the 4 compact images 
(images 1, 3, 4 and 6, filled squares). The model in the right panel is also constrained by the other 
images in the lens (the two-image system la/8, open squares; the four-image system 2a/2b/5/7 
filled triangles; and the two-image system comprising parts of 5/7, open triangles) The tangential 
critical line of the model (heavy solid curve) must pass between the merging images 2a/2b, but 
fails to do so in the first model (left panel). 
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Fig. 6. — Surface density contours for models of MG0414+0534 including misaligned as and 04 
multipoles (thin lines). The model in the left panel is constrained only by the 4 compact images 
(Ai, A2, B and C, filled squares). The model in the right panel is also constrained by the structure 
of the VLBI jets associated with each image. The heavy dashed curve is the tangential critical 
curve of the lens. Note the small satellite galaxy to the North of the lens - this is an example of 
the high mass substructures (M ~ IO^^Mq) which were included as part of the macro model in 
DK02 rather than treated as substructure. 
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3. Ruling Out Systematic Problems in the Model Potentials 

If the anomalous flux ratios cannot be explained by propagation effects in the ISM, the only 
alternative to some form of substructure in the gravitational field is a systematic problem in the 
assumed macro model for the lens galaxy. Previous studies of the effects of changes in the macro 
model show that changes in the radial structure from the standard singular isothermal ellipsoids 
(SIE) with external shear do not provide a solution (Metcalf and Zhao 2002, Keeton et al. 2002). 
This makes physical sense because the images in most of the systems with anomalous flux ratios 
all lie close to the Einstein radius of the lens, making it difficult for changes in the radial profile to 
produce changes in the flux ratios. We focus on changes in the angular structure because they can 
more easily shift the positions of the tangential critical line so as to change the image fluxes. 



3.1. Lens Models With Deviations from Ellipsoidal Structure 

We modified the gravitational potential by adding higher order multipoles of the form 



-r cos m{9 — Om) (1) 



m 

to the primary lens galaxy. These correspond to a surface density 

i^m = —- — ^ cos m(6' - (2) 
r m 

and match the expansion used by Evans &; Witt (2002). The new terms have been included in a 
revised version of the lensmodel package (Keeton 2001).^ In lensmodel, the surface density of the 
SIE model is defined by 

V2 r [1 + g2 _ (1 _ g2) cos 2{e - 62)^^ 

where g = 1 — 62 is the axis ratio and O2 and is orientation of the major axis. The critical line, 
which is an isodensity contour with n = 1/2 for the SIE, has a major axis of h and a minor axis 
of hq. Thus, if we add an m = 4 term to an SIE model, these normalizations imply a standard 
amplitude for the deviation of an isophote from the ellipse of = em(l ~ rn?')/mh. Positive £4 
corresponds to positive 04 and "disky" density contours when the m = 4 component is aligned 
with the quadrupolc (^4 = ^2)- Galaxies with "boxy" isophotes will have negative 04 and £4. 

We will restrict our attention to adding m = 3 and m = 4 terms to the potential of the 
primary lens galaxy. This should be sufficient to indicate their potential importance without 
the models being too under-constrained. Observational evidence indicates that the amplitude of 



^The lensmodel package now allows terms of the form {em/m)r" cos m{9 — Om) where the normalization was chosen 
to reduce to an external shear for n = 2 and m = 2. 
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these higher order multipoles is small. For example, Bender et al. (1989) and Rest et al. (2001) 

examine a sample of local ellipticals and quantify the deviations of the isophotes from ellipses. 
They find that the deviations are dominated by the m = 4 term, oc a4Cos(^ — ^4), with a typical 
amplitude of |a4| ~ 0.01. Deviations from ellipticity in projections of simulated halos are of 
comparable amplitude (Heyl et al. 1994, Burkert Sz Naab 2003). Considerably less is available on 
the distribution of misalignment angles, Om — O2, between the higher order multipoles and the 
dominant ellipse. We know from the statistics of lens models that the major axes of the mass 
and the light are closely aligned (Kochanek 2002). The effect of these higher order terms on 
lenses were first studied by Evans &: Witt (2001), who pointed out that larger values of 04 than 
are typically observed would lead to significant cross sections for lenses producing more than 4 
observable images. 

We fit the 7 lenses in the DK02 sample with a sequence of four models whose results are 
presented in Table 1. In the first model ( "astrometry" ) we fit only the image positions using our 
standard SIE plus external shear model. In the second model ( "astrometry+flux" ) we added the 
constraints from the image fluxes assuming 5% uncertainties. In either case, 6 of the 7 lenses have 
unacceptably high Xfix statistics for the image fluxes even though the standard models provide 
good fits to the astrometric data, Xtot ~ x'flx- Only B 1608+656 has a xjix largely consistent with 
5% measurement uncertainties. The failure of these models to reproduce the flux ratios is the 
origin of the entire problem of flux ratio anomalies. In the third model we added the m = 3 and 
m = 4 multipoles to the potential of the primary lens, but constrained them to be aligned with the 
ellipsoid (0„, = 62 but allowing any sign for e^). For MG0414+0534, B0712+472, PGl 115+080 
and B1933+503 the new terms lead to a significant reduction in Xjix^ ™ case docs it 
become a statistically acceptable fit. The amplitudes of the multipoles are unreasonably large 
for B0712+472 and B1933+503. Finally, in the fourth model, we allowed the m = 3 and m = 4 
multipoles to be misaligned with respect to the ellipsoid. In these models the flux ratio anomalies 
implied by the standard lens models for PG1115+080, B1422+231 and B1933+503 are gone, but 
the multipole amplitudes are still unreasonably large. 

Thus, like Evans &; Witt (2002), we find that adding higher order multipoles to the potential 
could explain some of the flux ratio anomalies (3 of the 6 in the DK02 sample) , if the positions and 
fluxes of the 4 compact images were the only available data. The surface density contours implied 
by these models are not very physical, as shown for PG1115+080 and B1933+503 in Figs. 4 and 
5. The improved fit with the higher order multipoles is easily understood for B1933+503. In this 
lens, image 4 is anomalously bright for an ellipsoidal model, and the new models solve the problem 
by making the tangential critical line lie much closer to image 4 so that it can be more highly 
magnified. The cause of the improvement in PG1115+080, where the anomaly is due to the flux 
ratio between images Ai and A2, is less obvious. In other cases, like MG0414+0534, the surface 
density contours, while very boxy, look more reasonable (Fig. 6). 

For three of the DK02 lenses with anomalies there are additional lensed components beyond 
the four compact images which can be used to test these models. In MG0414+0534 there are VLBI 
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jets associated with each of the quasar images (Ros et al. 2000, Trotter, Winn & Hewitt 2000), 
in PG1115+080 there is the Einstein ring image of the quasar host galaxy (Impey et al. 1998, 
Kochanek et al. 2001), and in B1933+503 there are additional multiply- imaged components of the 
radio source (Sykes et al. 1998, Cohn et al. 2001, Munoz, Kochanek & Keeton 2001). Table 2 shows 
the results of fitting the same series of models to these three lenses with the addition of the new 
constraints. For MG0414-I-0534, where the models were already unable to explain the anomalous 
flux ratios, little changes with the addition of the VLBI constraints, although the amplitudes of 
the multipolcs change as shown in Fig. 6. For the two lenses where it appeared possible for the 
higher order multipolcs to eliminate the anomalies, we find that the new constraints eliminate 
these solutions. Figs. 4 and 5 show the surface density contours for these new models, and it is 
immediately clear that the added constraints force the lenses to be significantly more regular and 
ellipsoidal. 



3.2. Statistical Checks of the Macro Model 

Thus, we find that higher order multipolcs can explain only 1 (B1422-I-231) of the 6 anomalous 
flux ratio lenses in the DK02 sample given the available data. Keeton et al. (2002) also noted that 
B1422+231 was less problematic than many of the other systems. Nonetheless, the perturbation 
amplitudes needed to explain B1422-I-231 are much larger than observed in real galaxies, halo 
simulations or in the lenses where additional constraints allow us to determine 03 and 04. We 
conclude that higher order multipolcs are not a likely explanation, but we would like a more 
generic test of this conclusion. The basic problem is that most lenses provide only enough data 
to constrain the standard ellipsoidal models, so more complex models that can explain the image 
fluxes are only constrained by the degree to which the models are viewed as physical. 

We noted in §2.2 that there is a statistical tendency for the anomalous flux ratios to appear 
as a suppression of the flux of the brightest saddle points relative to the predictions of smooth 
models, and that this was a characteristic property of low optical depth substructure. It is likely 
that any specific model which can explain the anomalies in individual lenses will be unable to 
reproduce this statistical property for randomly selected lenses. The worst case for examining 
this expectation is the effect of the m = 4 multipole, because it is the largest amplitude higher 
order multipole with a symmetry matching a common image configuration, the cruciform quads 
like Q2237-I-0305. If an m = 4 multipole is aligned with the ellipsoid it will preferentially affect 
the saddles as compared to the minima of a cruciform quad, but it will symmetrically affect both 
saddles and both minima. If it is randomly misaligned with respect to the ellipsoid, then it will 
not distinguish between saddles and minima, and if the observed lens geometries are dominated 
by merging pairs (fold catastrophes, like PG1115+080) or merging triples (cusp catastrophes, like 
B1422-I-231), then it will again be unable to distinguish saddles from minima. 

We tested these hypotheses for models consisting of an ellipsoid with an aligned (boxy or 
disky) a4 component or a randomly oriented 04 component for a range of models for the lens 
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Table 1. Basic Model Fits 



L(^i"is 


Case- 




- \ l 

A //.r 


A //.r 




In \\ 
1 " 1 1 


MG0414+0534 


astrometry 


2 


12.8 


41.2 








astrometry+flux 


5 


10.0 


35.5 








66^ = 


3 


2.4 


19.1 


0.011 


0.020 




full 


1 


6.3 


14.4 


0.005 


0.040 


B0712+472 


astrometry 


1 


4.2 


112.7 








astrometry+flux 


4 


31.8 


15.5 










2 


6.6 


21.4 


0.095 


0.178 




full 





1.9 


17.2 


0.242 


0.416 


PG1115+080 


astrometry 


2 


2.9 


57.0 








astrometry+flux 


5 


10.0 


40.3 








sem = o 


3 


2.2 


18.3 


0.006 


0.022 




full 


1 


2.9 


3.7 


0.022 


0.051 


B1422+231 


astrometry 


1 


2.1 


12.3 








astrometry+flux 


4 


2.7 


11.2 








sem = o 


2 


2.7 


11.3 


0.000 


0.000 




full 





0.5 


0.1 


0.089 


0.069 


B1608+656 


astrometry 


2 


1.8 


3.9 








astrometry+flux 


5 


1.6 


0.7 








50m = 


3 


2.2 


0.5 


0.000 


0.000 




full 


1 


4.1 


0.6 


0.000 


0.000 


B1933+503 


astrometry 


1 


0.5 


254.2 








astrometry+flux 


4 


19.2 


137.4 








sem = o 


2 


1.6 


0.1 


0.011 


0.137 




full 





1.3 


0.0 


0.079 


0.138 


B2045+265 


astrometry 


1 


1.4 


302.0 








astrometry+flux 


4 


2.1 


301.4 








sem = o 


2 


2.2 


301.4 


0.000 


0.000 




full 





1.9 


301.3 


0.001 


0.000 



Note. — Four model cases are fit: "astrometry" and "astrometry+flux" flt 
standard models (SIE + shear) to the image positions and fluxes. The "69jn = 0" 
adds the m = 3 and m = 4 multipoles constrained to be aligned with the ellipsoid 
{Om = 02), and the "full" models allow any orientation. We show the number of 
degrees of freedom Niiof, the goodness of fit Xtot ~ Xfix the constraints other 
than the flux ratios and the goodness of fit to the fluxes Xjix assuming 5% errors 
in the fluxes. The amplitudes of the multipoles are given by |a3| and |a4|. Very 
weak priors are included on the ellipticity and shear. 
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Table 2. Fits Testing the Basic Model 



Lens 


Case 


Ndof 


Xtot Xflx 


Xflx 


lasi 


04! 


MG0414+0534 


astrometry 


20 


52.7 


20.8 








astrometry+flux 


23 


53.0 


19.6 








66^ = 


21 


19.5 


17.6 


0.000 


0.032 




full 


19 


6.3 


15.2 


0.021 


0.027 


PG1115+080 


astrometry 


73 


222.1 


57.0 








astrometry + flux 


76 


222.7 


55.6 








dem = o 


74 


146.0 


13.8 


0.000 


0.027 




full 


72 


121.1 


22.9 


0.013 


0.021 


B1933+503 


astrometry 


11 


34.4 


268.2 








astrometry+flux 


14 


33.2 


256.1 








sem = o 


12 


44.7 


231.2 


0.007 


0.001 




full 


10 


33.5 


242.5 


0.009 


0.000 



Note. — Model sequences including higher order multipoles for the three 
lenses with flux ratio anomalies and additional model constraints. The columns 
are the same as in Table 1. 
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Table 3. Compton Limits on CLASS Source Sizes 



Lens 


S[m3y] 


u [GHz] 


z 


A*mod 


A6'inin[Mas] 


B0128+437 


18.9 


5 




4.3 


22 


MG0414+0534 


149 


5 


2.64 


18 


33 


B0712+472 


10.5 


5 


1.33 


17 


7 


B1422+231 


164 


8.4 


3.62 


8.2 


34 


B1555+375 


17 


5 




5.2 


18.6 


B1608+656 


34.1 


5 


1.39 


5.1 


23.8 


B1933+503 


17.6 


8.4 


2.62 


3.7 


14.7 


B2045+265 


29.02 


1.4 


1.28 


52 


24 



Note. — Minimum source sizes for radio lenses, based on the 
brightest image. Sources without redshifts were assumed to have 
z = 2. The Usted magnifications are from fits to the standard SIE 
+ shear model. We assumed a brightness temperature limit of 
Ti, = 0.5 X 10^2 
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population or typical lens geometry. If we simulate m = 4 models with the amplitudes typical of 
real elliptical galaxies or simulated CDM halos (|a4| ^ 0.02), the resulting ~ 5% flux perturbations 
are too small to be relevant for explaining the anomalous flux ratios. Wc must use models 
of abnormally large amplitude, |a4| = 0.05, to produce sufficiently large perturbations. After 
generating an image configuration, we fit the images using our standard SIE + shear models, 
fitting only the astrometric data and ignoring the fluxes. We then examined the statistics of the 
flux residuals after dividing the images into magniflcation and parity subsamples just as in §2.2. 
Wc focused on cusp and fold conflgurations dcflned by lenses with an opening angle between two 
images of < 30°. Of the DK02 sample, the systems MG0414+0534, B0712+472, PG1115+080, 
B1422+231, and B2045+265 would be classified as fold or cusp under this definition, while 
B1608+656 and B1933+503 just fafl to make the cut. 

The results of our Monte Carlo calculation are displayed in Figure 7a-c. As expected, 04 
perturbations to fold or cusp lenses do not systematically distinguish between image parities for 
disky, boxy, or randomly oriented m = 4 multipoles. This is in stark contrast to the properties of 
real fold and cusp lenses where the bright saddle points are preferentially demagnified compared to 
the other images (Fig. 7d). Clearly, m = 4 perturbations are incapable of producing this behavior, 
with KS test probabilities of less than 0.002 than any of these distributions agree with the 
distribution observed for the cusp/fold lenses in DK02. We experimented with still higher order 
multipoles (m > 4), and found that they also fare poorly in both systematically demagnifying 
saddles and differentiating between saddles and minima. The differences between distributions 
found in the Monte Carlo simulations and in the real data mean that the fiux ratio anomalies 
cannot be due to higher order multipoles. This would be true even if the models with higher order 
multipoles could successfully model the anomalous flux ratios in every lens! 



4. Ruling Out Microlensing of Radio Sources 

We are left only with substructure in the gravitational potential as a viable explanation of 
the phenomenon. Real lenses have substructures on two mass scales, stars and satellite halos, so 
the last step of our argument is to rule out stars as the source of the anomalous flux ratios of 
lensed radio sources. Microlensing is certainly an important phenomenon in the optical, where its 
effects are directly observed as uncorrelated time variability in Icnscd quasars (e.g. Q2237+0305, 
Wozniak ct al. 2000). The principal arguments against significant microlensing in radio lenses are 
laid out by Koopmans &: dc Bruyn (2000). The key point is that the typical scale of the caustic 
network is set by the Einstein radii of the stars {9e = 3[M /Mq]^/'^ j^as for a star of mass M, a lens 
at zi = 0.3 and a source at Zg = 2, hence the name microlensing), while the Compton limit for the 
minimum possible radio source size is given by 
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log(f,/f J 

Fig. 7. — Cumulative distributions of model flux residuals, log (/obs //mod)) expected from Monte 
Carlo simulations of fold and cusp lenses with m = 4 components of amplitude |a4| = 0.05. Panels 
(a) and (b) correspond to perturbations aligned with the major axis of the ellipsoid, for 04 > 
(disky) and 04 < (boxy) respectively. Panel (c) shows results for randomly oriented m = 4 terms. 
Panel (d) shows the flux anomalies observed for the 5 lenses in the DK02 sample satisfying the 
definition of fold/cusp described in the text. The curve labeling is the same as in Figures 2 and 3. 
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Table 3 shows that the minimum source sizes for the four-image CLASS lenses are typically 
^min ^ 10/ias, or roughly 10 times larger than for a stellar mass function dominated by low-mass 
stars. A sufficiently non-thermal electron distribution can change these limits, but it is unlikely 
to reduce them by an order of magnitude. While microlensing leads to fractional flux variations 
of order unity for sources smaller than the typical caustic size, larger sources average over the 
magnification pattern and show far smaller variations. For sources larger than 20/xas we would 
expect ^ 5% fluctuations from stellar microlensing. For microlensing to produce the ^ 20% flux 
anomalies of the radio lenses we would need to have source sizes of order ~ ^as. 

One means of escaping the Compton limit on the radio source sizes is to make most of the 
radio flux come from a relativistically beamed source. If a significant fraction of the source flux 
comes from a source moving towards us at relativistic speeds, then the brightness temperature of 
the beamed component can be increased by a Doppler factor V = [7(1 — /3cosi)]~^ ~ 7, allowing 
a corresponding reduction in the source size. For example, there are ~ 5% uncorrected changes 
in the radio fluxes of the lensed images in B1600-I-434 that can be explained by microlensing a 
beamed radio component with 7 10 (Koopmans & de Bruyn 2000) . Much larger Lorentz factors 
would be needed for beaming to explain the ^ 20% amplitude of the flux anomalies. 

In the absence of beaming, microlensing leads to changes in the image flux ratios on time 
scales of years, while CDM substructure would lead to changes only on time scales of millcnia 
or longer. If there is beaming, so that the characteristic velocities are superluminal rather than 
simply the peculiar velocities of the lens, source and observer, then the time scales will be 1000 
times shorter. Where radio sources have been monitored for extended periods (e.g. B1600-I-434, 
Koopmans & de Bruyn 2000; B1933-F503, Biggs et al. 2000; B1608-F656, Fassnacht et al. 2002) 
the average flux ratios show little time dependence. While this is only weak evidence against 
microlensing in the absence of beaming, it is strong evidence against microlensing if there is 
beaming because of the very short time scales created by the relativistic motions. Finally, if 
microlensing were important, we would expect a significant frequency dependence to the flux 
ratios because we expect the source size to be a function of wavelength. As discussed in §2.1, this 
frequency dependence has not been observed. 

Because of the stress we have put on the dependence of the flux residuals on image parity, we 
examined whether the differences between saddle points and minima we discussed earlier persists 
for source sizes comparable to the Einstein radii of the substructure. We ran inverse ray-tracing 
simulations (e.g. Schneider et al. 1992) using 10240^ rays and the particle-mesh (PM) force 
calculation algorithm (Hockney and Eastwood 1981). The PM method enjoys several advantages 
for the microlensing problem, chiefly in the simplicity and speed of the algorithm. Using this 
code, we simulated microlensing at a saddle-point and at a minimum, each with an unperturbed 
macro-model magnification of about 10. The saddle has a total convergence and shear (coming 
from the macro model) of kq = 0.525, 70 = 0.575, while the minimum has kq = 0.475, 70 = 0.425. 
We lay down a random distribution of microlenses with Einstein radii 6-^ = 1/ias and a number 
density corresponding to convergence = 0.075. The remaining surface density, kq — k^,, was in a 
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smooth component. 

In Figure 8 we plot the distribution of flux residuals as a function of source size. For a source 
small compared to ^e, like the case Vg = 0.13/xas corresponding to our grid resolution, we see that 
the magnification distribution for the saddle point is skewed towards demagnification compared to 
that for the minimum, in agreement with previous results (e.g. Schcchtcr & Wambsganss 2002; see 
their Fig. 3). For larger sources, which we illustrate for the cases = Oe = 1/uas and 39e = 3//as, 
the widths of the distribution narrow rapidly (with a^/iiQ cx l/r^ roughly, e.g. Koopmans &; de 
Bruyn 2000) and the differences between the distributions for saddle points and minima vanish. 
The mean (logio/^/Mo) provides a measure of the skewness, and for source sizes of Vg = 0.13, 
0.56, 1 and 3/xas we find values of (log^oM/z^o) = —0.26, —0.025 —0.01 and —0.002 respectively. 
Evidently, even source sizes as small as 9^/2 = 0.5;uas are sufficient to erase the asymmetry in 
the distribution of perturbations to saddle-point images. In reality, radio QSO's are expected to 
have much larger angular sizes (10-30/Ltas, see Table 3), giving even narrower, more symmetric 
distributions than those plotted in Figure 8. 

5. Summary and Discussion 

We have considered three alternatives to substructure as the source of the flux anomalies 
observed in four-image radio lenses: propagation effects like refractive scattering in the interstellar 

medium (ISM) of the lens galaxy, problems in the macro model describing the gravitational 
potential of the lens galaxy, and microlensing by ordinary stars in the lens galaxy. We find that 
none of these alternatives are likely to explain the phenomenon. 

Generating the flux anomalies using the ISM fails on two counts. First, the anomalies in 
the flux ratios of the radio lenses are essentially independent of wavelength, unlike any normal 
propagation effect in the ISM (scintillation, refractive scattering, free-free absorption). Any 
attempt to use the ISM as an explanation of the anomalous flux ratios must find a mechanism to 
produce a nearly frequency independent optical depth - the radio equivalent of the "gray dust" 
proposals for Type la supernovae (Aguirre 1998). Second, the fiux ratio anomalies depend on the 
image parity and magnification. In particular, the brightest saddle point image has a tendency 
to be demagnified relative to the other images, which is a characteristic of low optical depth 
substructure (Schechter &; Wambsganss 2002, Keeton 2002). Since the ISM should preferentially 
perturb the least magnified images, which will be the most compact, and show no dependence on 
the image parity, we can statistically rule out any propagation effect. 

Solving the flux anomalies by modifying the smooth gravitational potential of the lens also 
fails on two counts. Following other studies (Evans Sz Witt 2002, Moeller, Hewett & Blain 2002, 

Quadri, Moeller & Natarajan 2002), we investigated adding higher order multipoles {(j) oc cos mO 
with m > 2) to the potential. For the models we attempted (which were restricted to adding 
m = 3 and m = 4) we could find solutions without fiux anomalies for only one of the six lenses 
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log(/obs//mod) 

Fig. 8. — Cumulative distributions of magnification perturbations for microlensing, as a function 
of source size r^. The solid lines correspond to a minimum with kq = 0.475 and 70 = 0.425 while 
the dashed hnes correspond to a saddle point with kq = 0.525 and 70 = 0.575. The simulations 
used microlenses with Einstein radii of 6e = l^diS and surface density = 0.075. 



-23- 



from DK02 which had anomahes in our standard eUipsoidal lens models. The fits required the 

higher order multipoles to be significantly misaligned with respect to the axes of the ellipsoid and 
had perturbation amplitudes larger than seen in real galaxies or halo simulations. In three of 
the lenses, the presence of additional constraints, either Einstein ring images of the quasar host 
galaxy or additional lensed components of the radio source, allowed us to measure directly the 
magnitude of the higher multipoles in the gravitational potential. In every instance, the m = 3 
and 4 multipoles were constrained to the small amplitudes typically observed in galaxies and halo 
models, largely justifying the use of ellipsoidal mass models to describe lens galaxies. 

The second problem with the higher order multipole models is that they generally cannot 
reproduce the observed parity dependence of the flux anomalies. They can do so in very specific 
cases where the symmetry of the lens (a cruciform quad like Q2237+0305) matches the symmetry 
of the multipole (a cos 4^ multipole), and the multipole has very large amplitude and is aligned 
with the dominant ellipsoidal potential. However, none of the radio quads with anomalies have 
this symmetry (they generally have merging image pairs or triples associated with fold and cusp 
caustics) and the models require multipoles that are not aligned with the ellipsoid. If we examine 
the distribution of anomalies found by fitting lenses with higher order multipoles using only 
standard ellipsoidal lens models, we find that the statistical properties of the brightest saddle 
point do not stand out from those of the other images. This allows us to rule out these models as 
an explanation even if they could fit the anomalies in each individual lens (which they cannot). 

The failure of the ISM or changes in the smooth potential to explain the data means that 
the explanation must be substructure in the potential. These substructures must either be dark 
objects, since luminous satellites are not abundant enough to account for the preponderance of 
anomalous fluxes (e.g. Chiba 2002), or the stellar populations in the lens galaxy. Either source 
of substructure will lead to the statistical differences between saddle points and minima if the 
optical depth is low (Schechter Sz Wambsganss 2002, Keeton 2002). The angular scales of the 
magnification patterns are very different for stars (microarcseconds) and satellites (milliarcseconds) 
and by using radio sources we should be averaging out the stellar magnification patterns to leave 
only the contribution from the satellites. We tested this by taking typical microlensing patterns 
for saddle point and minimum images, which show the expected differences for point sources, and 
smoothing on larger and larger scales. As expected, the magnification distributions for the saddle 
point and the minimum show significant differences only when the source size is smaller than the 
Einstein radius of the stars. Source sizes of even one microarcsecond are sufficient to eliminate the 
differences, so radio sources could have angular sizes even an order of magnitude smaller than the 
Compton limit without microlensing making a significant contribution to the flux anomalies of 
radio lenses. 

In summary, CDM substructure remains the best explanation of the flux ratio anomalies of 
gravitational lenses. The most powerful piece of evidence is the statistically significant difference 
between the anomalies in saddle points and minima and between the bright and faint images. 
Substructure makes a very specific prediction that the brightest saddle point should show a very 
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different distribution from the other images, as observed in the data, that is very difficult or 
impossible to reproduce using either the ISM or changes in the smooth potential. 

However, the issues raised by these considerations point to future observations that can further 
clarify the origin of anomalous fluxes in lens systems. ISM effects can be more strongly constrained 
by measuring flux ratios at still higher frequencies (e.g. 43 GHz at the VLA). Mid-infrared 
(5-10/im) flux ratios, where the wavelength is far too short to be bothered by electrons and 
far too long to be bothered by dust, are difficult to measure but completely insensitive to the 
ISM (e.g. Agol, Jones &; Blaes 2000). The mid- infrared source sizes should also be large enough 
to be affected only by satellites, adding a further check for whether we can be misinterpreting 
microlensing effects as substructure. Integral field spectroscopy, to compare the flux ratios in 
the optical continuum to those in the broad emission lines, provides another way to separate 
the effects of microlensing and substructure because the broad line emitting regions should be 
significantly larger than the source of the optical continuum (e.g. Moustakas k. Metcalf 2002). 
Monitoring the lenses, either in the radio or the optical, to look for changes in the flux ratios, 
can also be used to distinguish microlensing effects from CDM substructure. Observations to find 
additional lensed structures are the best approach to determining whether more complicated lens 
potentials arc needed. Clean constraints on more complicated angular structures can be obtained 
by analyzing the shapes of the Einstein ring images of host galaxies found in deep, high-resolution 
infrared images (see Kochanek et al. 2001). 

We thank Roger Blandford, Scott Gaudi, Gil Holder, Chuck Keeton, Feryal Ozel and David 
Rusin for helpful discussions. We also thank Chuck Keeton for rapidly adding the higher order 
multipolc potentials to the lensmodel package, and Chris Fassnacht for giving us the unpublished 
5 and 15 GHz flux measurements for B1608-I-656. ND acknowledges the support of NASA through 
Hubble Fellowship grant #HST-HF-01148.01-A awarded by STScI, which is operated by AURA 
for NASA, under contract NAS 5-26555. CSK is supported by the Smithsonian Institution and 
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